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Polymorphism offers rich and virtually unexplored space for discovering novel functional materials. 
To harness this potential approaches capable of both exploring the space of polymorphs and assessing 
their realizability are needed. One such approach devised for partially ionic solids is presented. The 
structure prediction part is carried out by performing local DFT relaxations on a large set of random 
supperlattices (RSLs) with atoms distributed randomly over different planes in a way that favors 
cation-anion coordination. Applying the RSL sampling on MgO, ZnO and Sn02 reveals that the 
resulting probability of occurrence of a given structure offers a measure of its realizability explaining 
fully the experimentally observed, metastable polymorphs in these three systems. 


The discovery of polymorphism in the late 18th and 
early 19th century [1,2] revealed the significance of struc¬ 
tural degrees of freedom in determining physical proper¬ 
ties of solids. The best known example is probably el¬ 
emental carbon with markedly different mechanical, op¬ 
tical and electronic properties between its graphite and 
diamond forms [3]. Other notable cases include white 
and grey tin, which also exhibit significant differences 
in electronic and mechanical properties [3]; or enhanced 
photocatalytic activity of anatase Ti 02 compared to the 
ground state rutile polymorph [4]; or elemental Silicon, 
an indirect band-gap semiconductor in the ground state 
diamond structure predicted to become a direct gap ma¬ 
terial in a number of higher energy structures [5] includ¬ 
ing the experimentally realized clathrate structure [6] . 

However, the development of rational approaches to 
explore the space of polymorphs and (desirably) assist 
in their experimental realization faces significant chal¬ 
lenges. First, the complexity of the potential energy sur¬ 
face (PES) of periodic systems, evidenced by the expo¬ 
nential increase in the number of local minima with the 
system size [7] , limits our ability to systematically explore 
the spectrum of possible structures. 

A related problem of finding the ground state struc¬ 
ture attracted attention, especially with the development 
of first-principles total energy methods, resulting in a 
number of structure prediction techniques [8]. These 
include simulated annealing [9, 10], methods based on 
evolutionary algorithms [11-13], metadynamics [14, 15], 
basin and minima hopping [16, 17], random structure 
searching [18], methods based on data mining and ma¬ 
chine learning [19], structure prototyping [20, 21], etc. 
Although focused on finding the ground state structure 
some of these methods were also used in exploring the 
space of polymorphs (see for example Refs. [5, 22]) with 
the energy above the ground state as the main quantifier 
of their potential for experimental realization. 

This brings us to the second major challenge, the as¬ 
sessment of the likelihood for experimental realization of 
different polymorphs. While certainly being an impor¬ 
tant quantity, the energy above the ground state alone 



FIG. 1. (color online) A sketch of the potential energy sur¬ 
face (PES) of solids with different polymorphs corresponding 
to different PES local minima. 


is insufficient to explain observations based on available 
experimental data. Eor example, in the case of MgO, 
despite predictions [23] only the ground state rocksalt 
structure, and no other, is experimentally realized as re¬ 
ported in the Inorganic Crystal Structure Database [24] . 
In the case of ZnO, only the ground-state wurtzite and 
two other structures, zincblende and rocksalt, are exper¬ 
imentally realized [24-26]. Another important example 
is Sn 02 , which undergoes a series of phase transitions 
under pressure [27], but all of the high pressure phases 
relax to either the ground state rutile or the metastable 
(a-Pb 02 structure type upon releasing the pressure [28]. 
These facts indicate that for a given composition there 
seems to exist a finite set of structures that have higher 
likelihood for realization than the rest. 

If this is true, then the realizability of a given poly¬ 
morph can be thought of as determined by a combina¬ 
tion of three factors: (z) the energy above the ground 
state, {ii) the energy barrier to escape from a given PES 
minimum, and {in) the volume of configuration space oc¬ 
cupied by the PES minimum. The (z) and (zz) describe 
the principles of energy minimization and kinetic trap¬ 
ping. The factor (zzz) on the other hand, measures the 
probability of getting into a given PES minimum. More 
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FIG. 2. (color online) Steps in the random super lattice (RSL) structure generation (upper panel). Lower panel shows the 
number of distinct structures (polymorphs) within /\E = 50, 75,100,125,150 meV/atom above the ground state as a function 
of the total number of structures used in the sampling procedure. 


precisely, as shown in Fig. 1 every local minimum on the 
PES defines its basin of attraction, or the region of con¬ 
figuration space that has that minimum as its ’’center of 
gravity”. Hence, the probability of ’’falling” into a cer¬ 
tain structure has to be proportional to the total volume 
of configuration space occupied by its basin of attraction 
including all of the symmetry equivalent basins. 

Herein it is demonstrated, using MgO, ZnO and Sn02 
as case examples, that the factor {in) is actually critical 
in establishing a ranking of realizability with (z) and {ii) 
providing additional constraints. This is done by pursu¬ 
ing the idea that the total volume occupied by various 
basins of attraction can be measured using a large num¬ 
ber of random structures (random unit cell vectors and 
random atomic positions) that are relaxed to the closest 
PES local minimum utilizing density functional theory 
(DET). The frequency of occurrence of a given structure 
would then provide a measure of the probability to ” fall” 
into its local minimum. 

To do this and, at the same time, to overcome in part 
the difficulties posed by the already mentioned complex¬ 
ity of the PES, a structure prediction method is pro¬ 
posed to bias the random sampling toward the region 
of the PES more relevant for ionic systems. Because of 
the charge transfer only the structures that have cations 
preferentially coordinated by anions and vice-versa are 
relevant. The method adopted here favors the cation- 
anion coordination by distributing different types of ions 
in a random fashion over two interpenetrating grids of 
points. The grids are constructed using the alternating 
planes of a superlattice defined by a randomly chosen re¬ 
ciprocal lattice vector (see Eig. 2). Constructed in this 


way, these random superlattice (RSL) structures exhibit 
dominant cation-anion coordination. 

Eor each MgO, ZnO, and Sn02 a total of 2000 RSL 
structures with sizes varying between 1-20 formula units 
are constructed and DET-relaxed to the closest local min¬ 
imum. The relaxed structures are sorted into classes of 
equivalence and for the classes with largest occupancies 
(frequencies of occurrence) additional phonon calcula¬ 
tions are performed with the purpose of providing the 
information on the dynamic stability. The analysis of 
the resulting frequencies of occurrence shown in Eig. 3 
reveals that the experimentally observed polymorphs are 
exclusively the ones with the highest occurrence. 

Namely, for MgO the relative frequency of occurrence 
of the rocksalt phase is ^0.44, and is more than an order 
of magnitude larger than for any other structure, which 
explains why is the rocksalt phase so ubiquitous in MgO. 
Experimentally realized wurtzite, zincblende and rock- 
salt ZnO are the top occurring structures with the rel¬ 
ative frequencies in the 0.04 - 0.11 range well separated 
from the rest. In the case of Sn02, the relative frequen¬ 
cies are all below 0.02. The top two occurring are the 
ground-state rutile and the anatase structures. The a- 
Pb 02 is the structure type with the highest occurrence 
among those that have their total energy close to, and 
volume smaller than the rutile. These results elucidate 
why is it relatively easy to push Sn 02 outside of the rutile 
structure using pressure (relatively small basin of attrac¬ 
tion) and why it undergoes a series of phase transitions 
(all other structures have small basins), and why does 
(a-Pb 02 occur upon releasing the pressure. 

RSL Sampling. The details of the RSL structures gen- 
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FIG. 3. (color online) Relative frequencies of occurrence of structures resulting from the RSL sampling shown against the 
energy above the grounds state. Only the top occurring structures have their space groups explicitly marked. Space groups 
shown in grey represent the structures that are predicted to be the dynamically unstable. 


eration are shown in the upper panel of Fig. 2 . It is a 
modification of the Ab Initio Random Structure Search¬ 
ing method [18] and similarly starts with the random 
choice of unit cell parameters a, 6 , c, a,/3, 7 . In the sec¬ 
ond step the cation and anion grids are constructed in the 
following way. First, a transformation to the fractional 
(crystal) coordinates is performed to provide a cubic-like 
representation of the unit cell. Then, a reciprocal lattice 
vector G = nigi + n 2 g 2 + ^sgs with random ni,n 2 ,n 3 
is constructed. G defines a plane wave cos (Gr) and an 
associated superlattice. The two grids are constructed 
by discretizing the planes corresponding to the minima 
(cation grid) and the maxima (anion grid) of the plane 
wave. In the third step the ions are distributed over the 
two grids. To ensure homogeneous distribution and that 
that no two ions of the same kind are too close, the prob¬ 
ability distribution is constructed by placing a gaussian 
centered at each occupied grid point. The next ion is then 
placed on a grid point chosen randomly among those that 
have low probability. Finally, the structure is converted 
from fractional back into the real coordinates and the 
scaling factor is adjusted such that the minimal distance 
between any two atoms is larger than a certain thresh¬ 
old. The example of an RSL structure of ZnO shown in 
the Supplemental Material clearly displays the dominant 
cation-anion coordination. 

A total of 6000 RSL structures are generated (2000 
per system) with the following parameters: a, 6 , and c 
randomly chosen between 0.6 and 1.4 (in units of scale); 

and 7 random in the (30°, 160°) range; ni,n 2 , and 
ns also random between 4 and 10, the range which en¬ 
sures that sufficient, but not too large, number of planes 
in the unit cell; and the scale is adjusted such that the 
shortest distance between the atoms is not shorter than 
1.8 A. Different unit cell sizes are sampled by creating the 
RSL structures with one through 20 formula units and 
100 RSLs per size. Alternative would be to fix the cell 
size, but this would bias the sampling only to structures 
with sizes compatible with the chosen one. 

DFT Calculations. Full relaxations, including volume. 


cell shape and atomic positions, are performed on all 
RSL structures. This is done by employing standard 
DFT approach [29] with the PBE form of the exchange- 
correlation functional [30] and the projector augmented 
wave (PAW) method [31] as implemented in the VASP 
code [32]. The employed numerical setup (k-points, var¬ 
ious cutoffs) results in absolute total energies that are 
converged to within 3 meV/atom. All relaxations (vol¬ 
ume, shape, atomic positions) are conducted using the 
conjugate gradient algorithm [33]. Construction of the 
workflows, management of large number of calculations, 
and the analysis of results is carried out using the pylada 
software package [34]. 

Structure sorting. Sorting of the resulting, DFT- 
relaxed, structures into the classes of equivalence is done 
based on four criteria. Two structures are considered 
equivalent if: ( 1 ) their total energies are within 10 
meV/atom, (2) their space groups match, (3) their vol¬ 
umes per atom are within 0.5 %, and (4) the coordination 
of atoms up to the 4th neighbor is the same. After ex¬ 
tensive testing these four criteria were proven robust in 
establishing equivalence between the structures. 

Results. The RSL sampling procedure combined with 
DFT relaxations results in: 904 distinct structure types 
(classes of equivalence) for MgO, 1306 for ZnO and 1740 
for Sn 02 . The lower panel of Fig. 2 shows for all three 
systems the number of distinct structures within an en¬ 
ergy window AE = 50, 75,100,125, and 150 meV/atom 
above the ground state as a function of the total number 
of structures used in sampling. Two observations can be 
made. First, there are clear differences between the PESs 
of MgO, ZnO and Sn 02 , evidenced by a different num¬ 
ber of distinct structures present in the corresponding 
AE. Second, for all three systems the number of distinct 
structures grows monotonically with little, or no sign of 
convergence. This implies that even within a relatively 
narrow energy range a very large number of polymorphs 
can potentially exist. 

Eortunately however, a vast majority of the structures 
are actually irrelevant for their very low probability (fre- 
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quency) of occurrence. As shown in Fig. 3, except for 
a relatively small number of structures that occur more 
frequently, nearly all of the polymorphs resulting from 
the sampling procedure occur only once. Moreover, if 
the frequency of occurrence is represented as a function 
of the total number of structures, much like in the lower 
panel of Fig. 2, the extrapolation to the infinite num¬ 
ber of structures would imply zero probabilities for all 
of these (see Supplemental Material). All of the top oc¬ 
curring structures have their frequencies of occurrence 
fairly converged implying the finite probabilities of oc¬ 
currence. All of the experimentally realized polymorphs 
appear among the top occurring structures. 

In case of MgO, the top occurring is the rocksalt Fm- 
3m structure with the relative frequency of occurrence of 
0.437 far above any other structure. It is followed by the 
zincblende F-43m, and a four atom PGs/mmc structure 
with the relative frequencies of 0.0275 and 0.016, respec¬ 
tively. The occurrence of all other structures is below 
0.01. It is not a surprise then that the rocksalt is the only 
experimentally realized MgO structure. Note that these 
arguments do not mean that realizing other structures is 
entirely impossible, but that creating experimental con¬ 
ditions to do so might be much more challenging. 

The RSL sampling of ZnO polymorphs results in the 
distribution of the relative frequencies such that four dis¬ 
tinct structures occur more frequently than the others. 
These are: the ground state wurtzite PGsmc, zincblende 
P-43m, tetragonal /4mm, and the rocksalt Pm-3m. 
Phonon calculations reveal that the tetragonal /4mm 
shown in grey in Fig. 3 is dynamically unstable. The 
remaining three are exactly the three known, experimen¬ 
tally realized, polymorphs of ZnO [24-26]. The ZnO re¬ 
sults support directly the previous discussion about the 
energy above the ground state and its inadequacy to 
judge the realizability of different polymorphs. Namely, 
as shown in Fig. 3 there is a large number of distinct 
structures that appear inside the ^150 meV/atom en¬ 
ergy window between the ground state wurtzite and the 
rocksalt structure. A significant fraction of these have 
their volumes smaller than wurtzite. It is however, the 
rocksalt that is realized by applying the pressure despite 
its relatively high energy. Based on these results it can be 
argued that the realization of the rocksalt phase is due to 
the relatively large volume of configuration of space oc¬ 
cupied by the rocksalt structure. Of course, for the poly¬ 
morph to be metastable the kinetic barriers need to be 
sufficiently high to provide the conditions for the kinetic 
trapping, which is experimental fact for both zincblende 
and rocksalt ZnO [24-26]. However, it is the volume 
of configuration space that is comes first in determining 
the list of candidate structures with higher likelihood to 
be experimentally realized. The magnitude of the kinetic 
barriers will then determine which of these structures will 
actually be metastable. 

The results for Sn02 further support this discussion. 


Two most frequent structures in the RSL sampling are 
the ground state rutile P42mnm and the anatase lAiamd 
structure. As indicated in Fig. 3, cubic Fm-3m structure 
is dynamically unstable, while the fourth most frequent 
structure R-3m is ^200 meV/atom above rutile. The 
metastable Pbcn structure ((a-Pb02 type), that is ob¬ 
served in high pressure experiments, is the most frequent 
among the structures that have their energy close to, and 
volumes smaller than the rutile phase. This explains why 
is Pbcn Sn02 typically observed upon releasing the pres¬ 
sure. Interestingly, the anatase Sn02 predicted here to be 
among the most probable structures has so far not been 
realized experimentally. It has the volume significantly 
larger than the rutile and, contrary to Ti 02 which has the 
rutile and anatase polymorphs nearly degenerate, has the 
energy by more than 50 meV/atom above rutile Sn02. 
Therefore, the absence of anatase Sn02 is likely due to 
the challenges associated to achieving experimental con¬ 
ditions to fulfill both of these conditions simultaneously 
(expanded volume and higher energy). 

Finally, it is important to note that the actual values of 
the frequencies of occurrence do depend on the approach 
to sample different super cell sizes. As already mentioned, 
the supercells containing between 1-20 f.u. are sampled 
homogeneously in this work. The reason is to avoid bi¬ 
asing of the sampling only to sizes compatible with some 
fixed supercell size. Tests of different ways to sample the 
system size (fixed versus variable) show that the list of 
the top occurring structures that have their frequencies 
of occurrence converged with respect to the total number 
of structures does not depend on the sampling procedure. 

In conclusion, this paper demonstrates that the ran¬ 
dom superlattice (RSL) structure sampling followed by 
the DFT relaxations can be used to screen different poly¬ 
morphs of ionic systems as well as to assess the likelihood 
for their experimental realization. It is shown that the 
key quantity in assessing the likelihood for experimen¬ 
tal realization is the resulting frequency of occurrence, 
which measures (indirectly) the volume of configuration 
space occupied by a different structures. Application of 
the RSL sampling on MgO, ZnO, and Sn02 reveals that 
the experimentally observed polymorphs are exclusively 
the ones with the highest frequency (probability) of oc¬ 
currence explaining the physical reasons behind the ex¬ 
perimental observations. 
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